A condition number analysis of an algorithm for 
solving a system of polynomial equations with 
one degree of freedom* 

Gun Srijuntongsiri''^ Stephen A. Vavasis* 
December 21, 2009 

Abstract 

This article considers the problem of solving a system of n real polyno- 
mial equations in n + l variables. We propose an algorithm based on New- 
ton's method and subdivision for this problem. Our algorithm is intended 
only for nondegenerate cases, in which case the solution is a 1-dimensional 
curve. Our first main contribution is a definition of a condition number 
measuring reciprocal distance to degeneracy that can distinguish poor and 
well conditioned instances of this problem. (Degenerate problems would 
be infinitely ill conditioned in our framework.) Our second contribution, 
which is the main novelty of our algorithm, is an analysis showing that its 
running time is bounded in terms of the condition number of the problem 
instance as well as n and the polynomial degrees. 

1 Introduction 

We consider the problem of finding all zeros of a polynomial function / : 
[0, 1]"+^ M". The zero-set of such a function will, in the generic case, be 
a 1-dimensional algebraic set. Our algorithm is enumerative in nature, and 
therefore is feasible only in the case of small values of n. We refer to this 
problem as the single-degree-of- freedom polynomial system problem (SDPS). 

Perhaps the most common application of the SDPS problem is finding the 
intersection of two rational or polynomial surfaces, the so-called surface/surface 
intersection (SSI) problem. In this case, one is given two polynomials pi , p2 both 
mapping to K'^ that each parametrize a surfaces. The problem is to find their 
intersection, i.e., all points {s,t,u,v) in [0, 1]^ such that pi(s,t) ~ p2{u,v) — 0. 
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Other applications arise in robotics and motion planning. A final application is 
global optimization in which one finds all the local minimizers by constructing 
a network of paths that connect local minimizers (see, e.g., [T3]). 

Our proposed algorithm is a hybrid between subdivision and iterative meth- 
ods. This hybrid idea has been used to solve surface/surface intersection by 
Koparkar jT3] and to solve line/surface intersection by Toth [52]. The approach 
is to subdivide the domain recursively, discard the subdomains found to contain 
no solutions, and invoke an iterative method to locate a solution once it is cer- 
tain that the iterative method converges. The convergence tests used by Toth 
and Koparkar are both based on contraction mapping and evaluating ranges of 
functions. 

Our first main contribution, detailed in Section [3l is the definition of a 
condition number for SDPS problems. Intuitively, a problem instance is ill 
conditioned if it is close to a degenerate instance. A degenerate instance is one 
in which the Jacobian fails to have full rank at a root. Our condition number 
is the reciprocal of a quantity related to nearness to degeneracy. It is natural 
to expect that algorithms would have poorer behavior as the condition number 
grows larger. 

Our algorithm, which is presented in Sections is similar to Koparkar's in 
that it subdivides the parametric domains of the problem until the subdomains 
pass certain tests. It uses a bounding volume of a subdomain to exclude any 
that cannot have a solution. Our convergence test is based on the Kantorovich 
theorem, which tells us if Newton's method converges quadratically for the 
initial point in question in addition to whether it converges at all. For this 
reason, we can choose to hold off Newton's method until quadratic convergence 
is assured. Kantorovich's theorem is presented in Section [2] 

The main feature of our algorithm is that there is a lower bound on the 
size of the smallest hypercube occurring during the course of the algorithm that 
depends on the reciprocal condition number of the problem instance and on the 
polynomial degrees. Because our interest is in the low-dimensional and moder- 
ate degree case, we regard the factors depending on the degrees as 'constants' 
and the dependence on the condition number as the interesting feature. This 
analysis is presented in Section [T] A lower bound on the smallest hypercube 
size consequently implies an upper bound on the overall running time. 

As mentioned above, the SSI problem is a special case of SDPS with n = 3. 
Since there are many algorithms for SSI proposed in literature, we discuss here 
the main advantage of our algorithm as compared to other SSI algorithms. 
To the best of our knowledge, there is no previous algorithm for SSI in this 
class whose running time has been bounded in terms of the condition of the 
underlying problem instance, and we are not sure whether such an analysis is 
possible for previous algorithms. Indeed, we do not know of any SSI algorithm 
in the literature that has any a priori bound on the running time. We do know, 
however, that some algorithms do not have this property — their running time 
can be arbitrarily large even if the input instance is well conditioned. This is 
because these algorithms can sometimes create degenerate or nearly degenerate 
subproblems even though the original input instance is well conditioned. Section 
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[3] shows in details how marching methods based on coUinear normal points and 
Koparkar's algorithm in particular have the capacity to create bad subproblems 
from a good instance. 

The notion of bounding the running time of an iterative method in terms of 
the condition number of the instance is an old one, with the most notable exam- 
ple being the condition-number bound of conjugate gradient (see Chapter 10 of 
|10| ) . This approach has also been used in interior-point methods for linear pro- 
gramming [3], Krylov-space eigenvalue computation [5T], and the line/surface 
intersection problem [12] . We note that a related problem of computing convex 
hull of points on a plane is shown to always be well-conditioned [TT] . 

An additional motivation, not pursued further herein, for defining a condi- 
tion number and condition-aware algorithms like ours is that this creates the 
possibility of preconditioning. Preconditioning, which has been very successfully 
applied in numerical linear algebra (see, e.g., [23), means improving the con- 
dition number of an instance via some kind of transformation prior to solving 
it. 

We now define the problem under consideration more precisely by specify- 
ing a representation for the input polynomial system. Let Zi^rn{t) denote the 
Bernstein polynomials 



We are interested in finding all points x = {xi, X2, ■ ■ ■ , Xn+i) G [0, l]"^'^ satis- 
fying 



where ,i„_^i G M" {ij = 0, 1, . . . , nij) denote the coefficients, also known as 
the control points. Therefore, the problem is specified by the (mi -I- l)(m2 -I- 
1) • • • (m„-|_i -|- 1) control points. (See a further remark on this matter in Sec- 
tion [5]). This form of a multivariate polynomial is sometimes called tensor 
product Bezier representation: it presumes that the maximum degree of vari- 
able Xi (separately) is for each i = l,...,n-|-l. Note that / is a function 
that maps R^^^ to M". Note that this representation would be intractable for 
a large value of n, but, as mentioned earlier, our algorithm is intended for small 
values such as n = 3. 

The Bernstein basis is known to have better numerical stability for poly- 
nomials on the unit interval than the power basis [HI [7] , and computation us- 
ing parametric representation is often much more efficient than other types of 
surface representations. Furthermore, our algorithm makes direct use of the 
Bernstein-Bezier representation. In particular, our exclusion test is based on 
Bezier control points. It should be noted that the algorithm proposed in this 
article can be generalized to use with parametric surfaces represented by other 
polynomial bases provided that an appropriate exclusion test is available, and a 
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few other properties hold for the basis. Refer to [in] for a related algorithm for 
line/surface intersection that can operate on parametric surfaces represented by 
other polynomial bases. 



2 The theorem of Kantorovich 

Denote the closed ball centered at x with radius r > by 

B{x,r)^{yeW' : <r}, 

and let B{x, r) denote the interior of B(x, r). Kantorovich's theorem in aflinely 
invariant form, which is valid for any norm, is as follows. 

Theorem 2.1 (Kantorovich, affinely invariant form [51[T^). Let f : D C M" 

M" be differentiable in the open convex set D. Assume that for some point 
G D, the Jacobian is invertible with 

\\f{x°)-'f{x^)\\<V- 
Let there be a Lipschitz constant u > for f {x'^)^^ f such that 

<c..||x-y|| forallx,yeD. 
If h = r]LU < 1/2 and B{x°,p^) C D, where 

1 - Vi - 

P- = , 

Ll) 

then f has a zero x* in I3[x^ , p-). Moreover, this zero is the unique zero of f 
in {B{x°,p-)U B{x°,p+))r\D where 



1 + Vl-2/i 
P+ = 

and the Newton iterates x^ defined by 

x^+^ = x'' - fix'^y^fix'') 

are well-defined, remain in B{x'^ , p^), and converge to x* . In addition. 

We call x^ a fast starting point if the quantity h defined above satisfies 
h<l/4 and B(a;",p_) C D. In this case, quadratic convergence of the iterates 
starting from x*^ is implied. 
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3 A condition number of a polynomial system 
with one degree of freedom 



In this section we propose a definition for a condition number of tlie SDPS 
problem and prove tliat our condition number is related to tire distance from 
degeneracy. Let M denote the maximum p-norm among the control points of 
/ for some p. Later on, we will specialize to the infinity norm. It is easy to 
show that this quantity satisfies the axioms of a norm, so we will write M also 
as ||/||. Define the condition number oi f to be 

cond(f) = il/- max | min i - — ^^7-A\f'(x) 

.£[0,1]"+^ V iii/(a;)ir" 

Here, the notation means = (AA^)"^ . In the case that the rank of 
A is n, this corresponds to the Moore-Penrose pseudo-inverse of A, typically 
denoted as A'^. In general, the Moore-Penrose pseudo- inverse is defined for 
matrices of all ranks. In this paper, however, we need A^ only in the case 
that rank(()A) — n; we will take the second factor || of ([2|) to be oo 

if the rank of f'{x) is less than n. Similarly, the first factor is taken to be 
oo if f{x) = 0. Note that small condition number means the problem is well- 
conditioned. Indeed, it follows from ([3]) and (|16p below that both terms in the 
min have a lower bound of const /M. 

The rationale for this definition is that, as mentioned earlier, a degenerate 
instance has a point x such that f{x) = and rank(()/'(a;)) < n. For such 
a point, both terms occurring in the min of ([2]) are infinity, i.e., the condition 
number is infinite. Thus, the problem is ill-conditioned if f{x) is close to zero 
and f'{x) is close to rank-deficiency at the same point x. 

The inclusion of the factor of M makes the definition scale-invariant. Note 
that M depends on the choice of basis, namely, tensor-product Bernstein-Bezier 
basis, whereas the other factors in the condition number are basis-independent. 
The dependence on basis, however, is only up to a scalar factor depending on 
degree. This is because polynomials are a finite-dimensional vector space, hence 
all norms are equivalent. We could redefine M in a basis- independent manner 
as follows: 

Mbi= max 

xe[o,i]"+^ 

It follows from ^ and ([T7)) below that Mbi and M differ by a scalar that is 
bounded in terms of the degrees. The definition Mbi, however, would be harder 
to compute in practice. 

This condition number is similar to the definition of of Cucker et al. 
[3] for the problem of computing isolated roots of polynomial systems M" 
M". Our definition is somewhat simpler than theirs, however, because of the 
assumption we have imposed that the domain of interest is [0, l]""*"^ rather than 
all of R""*"^. This simplifying assumption obviates the need for introducing 
projective space and scaling of coefficients as in [3] . 
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The classical Turing Theorem [2] states that the condition number of a 
matrix, which bounds the iteration count of the conjugate gradient algorithm, 
is exactly the reciprocal of the relative distance of the matrix to singularity. 
Similarly, Shub and Smale show that their condition number for a homogeneous 
polynomial system is equal to the distance of the system to singularity [T7] . We 
now derive a result showing that our proposed condition number is also related 
to the distance of an SDPS instance to degeneracy. 

Before stating and proving the theorem, we require two well known bounds 
concerning polynomials in Bernstein-Bezier form: 

11/(^)11 < 11/11 (3) 

for all X € [0, 1]"+^. This follows because every value of / over the parametric 
domain is a convex combination of control points [6]. Next, 

||/'(a;)||oo < ll/'lloo < 2(n+l)max(mi,m2,...,m„+i)||/||oo (4) 

for all a; € [0,1]"'^^. This follows because the control points of a column of 
/' (i.e., a partial derivative of /) are finite differences of control points of / 
multiplied by the degree in the direction of differentiation as in (fT^ . A factor 
of 2 comes from the taking of finite differences, and a further factor of n + 1 
arises from the fact that the infinity norm is a sum over rows (not columns) of 
the derivative. Furthermore, applying the equivalence of norms to Q yields 

ll/'ll < 2Q:„(n + l)max(mi,m2,...,m„+i)||/|| (5) 

for any arbitrary p-norm, where a„ > is a scalar depending on n and the 
choice of norm. 

Now, finally, we come to the main theorem of this section. 

Theorem 3.1. Let f : [0, 1]"+^ — > M" be a polynomial function of degrees m = 
(toi, m2, . . . , rrin+i) in its n + I variables, and assume that it is nondegenerate, 
i.e., there is no x Q [0,1]"+"'^ such that f{x) = and rank(()/'(x)) < n. Then 
any f satisfying 

II / ~ /I ^ '^m /„\ 

11/11 - cond(/) ^ ' 

is nondegenerate, where Cm is a scalar depending on the degrees and the choice 
of norm. 

Conversely, there exists a degenerate polynomial f such that 

11/ /I ^ 

11/11 -cond(/)' 
where is another scalar depending on the degrees. 

Remark 1. As mentioned above, the vector and matrix norms appearing in 
([U, ([7]) may be any of the standard p-norms, although later we will specialize to 
the infinity norm. Inequalities ([6]) and ([7|) involve the norm of the polynomial 
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function. As mentioned above, we take this to mean the maximum norm control 
point when written in Bernstein-Bezier form. 

Remark 2. Let T) denote the set of degenerate polynomials, that is, those 
polynomials / such that cond(/) = oo. This theorem shows that our condition 
number is, up to constant factors, the reciprocal distance of / to 2? scaled by the 
norm of /. Consider a larger set V of degenerate polynomials defined as follows. 
Polynomial / g P' if there is any point x in C"+^ (not merely [0, 1]"+^) such 
that f{x) = and f'{x) is rank deficient. This set V is an algebraic variety, 
i.e., the coefficients of such /'s are the roots of a polynomial system. This 
means that we can apply Demmel's theorem [3] to conclude that the expected 
logarithm of the condition number of a random instance is modest. (Clearly, 
the distance of a polynomial / to I? is bounded below by the distance of / to 
V .) We can also apply the more recent analysis of Biirgisser et al. [1] to show 
that the 'smoothed' condition number [18j is modest, i.e., for any polynomial 
/ (even a degenerate one), if we select a random small perturbation of it, then 
the resulting polynomial is expected to have a modest logarithmic condition 
number. 

Proof. Let / satisfy (O and let e = / — / (i.e., a polynomial), so that ||e|| < 
Cm II /II / cond(/). Choose an x e [0, 1]"+^. By definition of cond(/), 

Af •min(l/||/(x)||,||/'(a;)t||)<cond(/). 

We take two cases depending on which term achieves the min. First, sup- 
pose M/||/(a;)|| < cond(/). In this case, ||/(a;) - /(a;)|l = ||e(a;)|| < ||e|| < 
Cm||/||/ cond(/). Assume Cm is sufficiently small so that < 1/3. Then 



< 



1 



1/(^)11 " \\f{x)\\-\\f{x)-f{x)\\ 

1 

< 



|/||/cond(/)-c,„||/||/cond(/) 
< 1.5cond(/)/||/||. 

In particular, f{x) ^ 0. 

For the other case, the hypothesis is ||/|| • ||/'(a;)t|| < cond(/). We recall 
that 



^ <||/'(x)t||o,<^;^±l (8) 



a„(/'(x))V^ - w noo _ 

(see [lOj (2.3.11) and §5.5.4]), where cr„(-) is notation for the nth singular value 
of a matrix. By the equivalence of norms, ([5]) is equivalent to 



7=<II/Wll<""^ 



for any arbitrary p-norm, where 5„ and 7„ are positive constants that depend 
on n and the choice of norm. 
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A second fact from numerical linear algebra is that 

\a,{A) ~a,{B)\<\\A- <V^\\A- B\\ 

(see [ini (2.3.11) and Cor. 8.6.2]). Combining these facts together with ([5]) yields 

C7,,{f\x)-e\x)) 
<^n{f\x)) - V^lle'll 

(Tn{f'{x)) - 2anVn{n + l)max(TOi, . . . ,m„+i)||e|l 

Sn 2any/n{n + 1) max(mi, . . . , mn+i)cm||/|| 

" y^|i/'(x)t|| cond(/) 
> ^«H/I _ ^an^/njn + 1) max(mi, . . . , TO„+i)cm||/|| 
Vf^cond(/) cond(/) 

If we select Cm < i5„/ (2Q!„n(n + 1) max(mi, . . . , m„^i)), then we are assured in 
this case that cT„(/'(a;)) > 0, i.e., the rank of /'(x) is n. 

Thus, combining the cases, we have shown that for all x G [0, 1]"^^, either 
/(x) 7^ or the rank of f'{x) is at least n, thus proving that / is not degenerate. 

Next, let us turn to ([7]). In this proof, we will drop the prime from and 
indeed will allow Cm to denote a constant depending on the degrees that may 
change from line to line. 

Let xq € [0, 1]"+^ be the point where the max in ^ is achieved, and let 
k be the value of this max, i.e., the value of cond(/). This means, first, that 
||/(a:^o)|| < M/k. Second, it means that either (i) |l/'(a;o)'l'|| > k/M, or (ii) 
rank(()/'(xo)) < n. Let A = f'{xo). 

If case (i) is true, then by definition of the matrix p-norm there exists a unit 
vector (in the norm under consideration) r such that || A^(ylA^)~^r|| > k/M. 
Let q = A^iAA^y^r, so that ||q|| > k/M and Aq = r. Let S ^ rq^ /{q^q) (a 
n X (n + 1) matrix) and define f{x) = f{x) — f{xo) — S{x — xq). In case (ii) 
(when rank(()/'(a;o)) < n), define f{x) = f{x) — f{xo), i.e., take 5 = 0, and we 
do not need r and q. 

For either case, we now must establish that / is a degenerate instance and 
that inequality ([7]) holds. 

First, let us establish the degeneracy of /. Clearly xq is a root of /. In 
case (ii), /'(xo) = /'(xq), a matrix whose rank is less than ri. In case (i), 
/'(xo) = /'(xo) — S = A — rq^ /{q^q). We claim that /'(xq) has at least two 
independent vectors in its null space, which implies that its rank is at most n—\. 
Observe first that A, as an n x (n + 1) matrix, must have a nonzero vector w 
in its null space. Then w is also in the null space of /'(xq) since Aw — and 
q^w — {AA^)^^ Aw = 0. Also, /'(xq) has q in its null space as f'{xo)q 
evaluates to Aq — r. Finally, q and w are independent since Aw = whereas 
Aq = r, which is not zero. This concludes the argument that / is degenerate. 

Next, consider \\f — f\\ appearing in ([7|). Let us introduce the following 
notation for this argument: r denotes the polynomial x i-^ x — xq and k denotes 
the constant real- valued polynomial 1. With these definitions, f — f = /(xo)k-I- 



^n(/'(a:)) = 
> 
> 
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St. Note that ||t|| < Cm since |jxo|| < 1. Also, = 1. We have the following 
chain of inequalities, which applies to case (i). The first line involves norms of 
polynomials, whereas the remaining lines are norms of vectors and matrices. 



< 
< 

< 

< 



||/(xo)k + S ot 

Cm (||/(a;o) 
'M 



< 



< 



M 

~k 

M 

T 

M 

~k 



+ 11^11) 
rq^/{q^q)\\ 

Hl-i/lkl 

r|| • M/k 



M 



(9) 



The first line follows from the definition of /. The second uses the fact that 
||k|| and ||r|j are bounded by constants. The third line uses the inequality 
established earlier that ||/(xo)|| < A4/k and also expands the definition of S. 
The fifth line uses the fact that q^q ~ WqW'ii a-nd the 2-norm and the p-norm 
under consideration are related by constants depending on mi, . . . , m„+i. The 
last line uses the assumption that r is a unit vector. 

For case (ii), (O also holds since 5' = 0, so the result is already established 
by the third line in the above chain of inequalities. Since the denominator 
occurring in the left-hand side of ([7]) is M, and recalling that k = cond(/), we 
see that ^ proves the theorem. □ 

Although the condition number defined by ^ is scale invariant (i.e., 
cond(/) = cond(c/) for c 7^ 0), it is not affinely invariant. In other words, 
if A is a nonsingular n x n matrix, then in general cond(/) 7^ cond(A/). On 
the other hand, our algorithm is affinely invariant as we shall see in Section [5l 
Therefore, we can define a new condition number that is indeed affinely invari- 
ant, which is as follows: 



icond(f) inf cond(Af). 

^ ' AeGL{nM) ' 



(10) 



Here, GL(n,R) denotes the set of all n x n nonsingular matrices. Obviously, 
pUj) is affinely invariant, and also, it is clear that for any instance / of SDPS, 
icond(/) < cond(/). Furthermore, if we are able to show that our affinely 
invariant algorithm has running time bounded in terms of cond(/), then it will 
follow automatically that it is also bounded in terms of icond(/). 

The difficulty with (fTO)) is that there is no obvious way to compute this 
quantity other than the exhaustive method of trying out all choices of A in a 
dense grid lying in GL(n,R). (If the matrix 2-norm is used, then it suffices to 
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try a dense sampling of upper triangular matrices, a smaller search space, since 
the Q in a QR factorization of A does not affect the norm.) Unless a better 
method can be found, definition (jlOp would be useful in practice mainly in cases 
where there is a priori information about a linear transformation that improves 
the condition number. 

4 Performance of other algorithms for n = 3 case 
on well conditioned instances 

Recall that SSI is a special case of SDPS with n = 3. Due to the abundance 
of SSI algorithms in literature, we compare our algorithm to well-known SSI 
algorithms. Many previously published SSI algorithms work well in practice 
and are widely used in computer-aided geometric design software. Nonetheless, 
we suspect that most of these algorithms can behave nonrobustly in the sense 
that, given a well conditioned problem instance, they can sometimes internally 
generate an arbitrarily ill conditioned subproblem that they then must solve. 
If an algorithm is capable of this behavior, then it is not possible to bound 
its running time in terms of the instance's condition number as we shall do 
for our algorithm. Indeed, as far as we know, there is no a priori running time 
upper bound for any SSI algorithm in the literature. In this section, we consider 
how two well-known SSI algorithms can generate bad subproblems given good 
problem instances. 

4.1 Marching methods based on colhnear normal points 

Collinear normal points are points on the two surfaces whose normals are 
collinear. Marching methods based on collinear normal points split the para- 
metric domains in at least one direction at these collinear normal points. The 
consequence is all solution curves have one point on the boundaries of the re- 
sulting subdomains provided that the dot product of any two normal vectors 
of either surfaces is never zero. These points are located by a curve/surface 
intersection algorithm and used as starting points for the marching step. 

Consider applying collinear normal points marching methods to find inter- 
sections between the two Bezier surfaces p and q whose control points are defined 
in Table[T]and Tabled respectively. This problem is equivalent to solving SDPS 
with control points 

where ai^^ij's (ij — 0, . . . ,mj) are control points of p and a^g^i^'s (ij — 0, . . . ,mj) 
are control points of q. The condition number of this instance is 423.4, which 
is reasonably well-conditioned. The instance is also intuitively well-conditioned 
as there are neither complicated nor almost singular intersections. The two 
surfaces and their intersection in object space are shown in Figure [TJ 

Figure [5] shows a pair of collinear normal points in parametric space and 
the subdivision of domain at these points, as well as the intersection between 
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(.155, .055, .002524)'^' 
(.655, .055, .005414)^ 
(1.155, .055, -.01745)^ 


(.155, .555, .003592)' 
(.655, .555, -.01454)^ 
(1.155, .555, -.02108)^ 


(.155, 1.055, -.008142)'^' 
(.655, 1.055, .005146)^ 
(1.155, 1.055, .01718)^ 



Table 1: The control points of Bezier surfaces p. The entry in the iith row and 
the «2th column is the control point a^j^ij of p. 
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2 
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(.1768, .1295, .02303)' 
(.4081, .1287, .04006)^ 
(.7515, .1249, -.07777)^ 
(1.6068, .1078, -.004634)^ 


(.1767, .3465, .06306)'' 
(.4081, .3434, .0948)^ 
(.7515, .3294, -.1071)^ 
(1.6068, .2651, .01833)^ 


(.1767, .6467, -.0801)' 
(.4081, .6384, -.08438)^ 
(.7515, .6008, .05544)^ 
(1.6068, .4288, -.01627)^ 


(.1768, 1.437, -.01946)'' 
(.4081, 1.444, .004442)^ 
(.7515, 1.477, -.01756)^ 
(1.6069, 1.628, .2468)^ 



Table 2: The control points of Bezier surfaces q. The entry in the iath row and 
the i4th column is the control point of q. 




Figure 1: An example of a well-conditioned instance of SSI problem where 
collinear normal point methods and Koparkar's algorithm require long compu- 
tation time to solve. 
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Figure 2: Intersection of the surfaces in Figure [T] in parametric space and split- 
ting of the domains at a pair of coUinear normal points showing that marching 
methods based on collinear normal points create ill-conditioned subproblem 
instances from a well-conditioned original instance. Solid lines are the intersec- 
tions. Asterisks are the pair of collinear normal points. Dotted lines are the 
splitting lines at those points. To the left is the parametric domain of p. To the 
right is that of q. 

the surfaces (Other pairs of collinear normal points are not shown). Parts of 
the intersection curves lie almost parallel to the boundaries of the created sub- 
domains. In other words, the cuts made by the algorithm happen to intersect 
degenerately or nearly degenerately with the problem data. Thus, after the al- 
gorithm has made a subdivision of the domain based on the collinear normals, it 
must now recursively solve arbitrarily ill conditioned subproblems since the in- 
tersection point between the surface and the boundary curve is almost singular. 
The running time of algorithms for finding these intersection points typically 
depends on their conditioning (see, e.g., the RIA algorithm of [E]). Thus, the 
running time of the algorithm cannot be bounded in terms of the condition num- 
ber of the input. Furthermore, ill-conditioning of the subproblems may cause 
unexpected inaccuracy of the solution of the original well-conditioned instance. 

4.2 Koparkar's algorithm 

Koparkar's algorithm uses a test based on the contraction mapping theorem 
to determine if, in a given domain, a Newton-like method converges or the 
two surfaces do not intersect at all. If the Newton-like method is guaranteed 
to converge, the method is used to find part of the solution curves inside the 
domain. If the two surfaces are known not to intersect, the domain is discarded. 
Otherwise, each surface is subdivided by splitting their parametric domains into 
four, and the test is repeated on the created subdomains. This process continues 
until the entire domain is examined. 
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Figure 3: Intersection of the surfaces in Figure [T] in parametric space and sub- 
domains created after the original domain [0, l]** fails Koparkar's test. Due to 
part of the intersections lying on a boundary of a subdomain, Koparkar's algo- 
rithm requires many subdivisions before those parts of the intersections can be 
located. Solid lines are the intersections. Dashed lines are the subdivision lines. 
To the left is the parametric domain of p. To the right is that of q. 

The test in Koparkar's algorithm requires the ability to evaluate ranges of 
functions, which is typically accomplished by variety of interval arithmetics. 
None of these techniques can give the exact ranges, however; they yield super- 
sets of the ranges. For this reason, the convergence test is likely to fail when 
part of the solution lies very close or directly on the border of a subcube in both 
xiX2-space and a;3X4-space at the same time, which is not necessarily on or near 
the border of the original domain [0, 1]"*. The same problem instance discussed 
in Section 14.11 which is shown in Figure [1] has such problem. The domain 
[0, l]** does not pass the test, and the domain is subdivided at the midpoints as 
shown in Figure [31 The subdomains now have solutions directly on a bound- 
ary. Koparkar's algorithm needs to subdivide these subdomains to very small 
ones before the convergence test is satisfied. This example demonstrates that 
Koparkar's algorithm is inefficient at solving certain well-conditioned instances. 
The algorithm may solve other instances with higher condition numbers but 
without any parts of the intersections near any boundaries of the subdivided 
subdomains faster than this well-conditioned instance. 
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5 The Kantorovich-Test Subdivision algorithm 



This section describes our algorithm for the SDPS problem. Some details are 
postponed to the next section. Since we are interested in solutions of / within 
the hypercube [0, 1]"'*'^, and the closed ball B{x, r) defined in the infinity norm 
is a hypercube, our algorithm uses the infinity norm for all of its norm compu- 
tation. Therefore, for the rest of this article, the notation ||-|| is used to refer 
specifically to infinity norm. 

During the computation, our algorithm maintains a list of explored regions 
defined as parts of the domain [0, 1]"+^ guaranteed by Kantorovich's Theorem to 
contain only the solutions that have already been found. This list is used in ad- 
dition to another test to determine whether to subdivide a hypercube. We define 
the Kantorovich test on a hypercube X = B{x'^,r) as the application of Kan- 
torovich's Theorem on the point x^ to the function h'^^^'>{x) — {f{x),Xi — k) 
for each i — l,2,...,?i + l and any k E [x^— r, .x^ + r]. The hypercube 
.5, 1.5]"+^ is used as the domain D in the statement of the theorem, and 

(a;0) h<-'''\x°) is used as 77. For lu, we instead use uj > uj, where 

a) is defined by (jl3p below, as the minimal lu is too expensive to compute. The 
hypercube X passes the Kantorovich test if there exists an i G {1, 2, . . . , n -I- 1} 
such that for every k E [x^ — r, x^^ + r], rjuj < 1/4 and B{x'^, p_) C D. 

The choice of D mentioned in the previous paragraph is used in the analysis 
of the algorithm in Section [71 but in practice a smaller D such that B{x''\ r) C 
D C [—.5, 1.5]"+^ may be used. The advantage of using a smaller D is that the 
Lipschitz constant d) will be smaller, so the inequality rjili < 1/4 may be satisfied 
more easily (i.e., for larger hypercubes) than using the full D. The disadvantage, 
however, is that if D is chosen too small, then the condition B{x^,p^) C D oi 
the theorem may be hard to satisfy. In particular, the choice D — B(xP,r) 
would be unacceptable for this reason. 

If X passes the Kantorovich test, then three important consequences follow. 
First, x° is a fast starting point for /i^*'^) for the particular i that satisfies the 
condition of the Kantorovich test and any fc € [x^ — r,x^ + r]. Second, the 
segment of the solution curve of / that contains the root guaranteed by the 
conclusion of Kantorovich's theorem is not a loop in 2:1X2 • ■ ■ x„-f i-space inside 
X (although it may be part of a loop in the original domain). Third, an explored 
region for this segment of the solution curve can be derived. The explored region 
is 

XE = {x:x'^-r<x,<x'^+r}nDn f] (b(xO, pL'' ^ S(xO, pf^^^)) , 

ke[x°-r,x°+r] 

(11) 

where p^'' and p^^'' are p_ and p_|_ in the statement of Kantorovich's theorem 
with respect to ft.'^'*'-' . Observe that Xe is a hyperrectangle in R"+^ and can be 
stored and computed succinctly as detailed in Section 16.21 Note also that the 
explored region provides an effective way to prevent the points on different but 
nearby solution curves from being incorrectly joined into the same curve. 
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The other test our algorithm uses is the exclusion test. For a given hyper- 
cube X, let fx be the Bernstein polynomial that reparametrizes with [0, 1]"+^ 
the function defined by / over X. In other words, fx{q) = where X{q) 

is a composition of a dilation and translation (uniquely determined) such that 
A : [0, 1]"+^ — + X is bijective. (See Section [6731 for information how to efficiently 
compute fx.) The hypcrcube X passes the exclusion test if the convex hull of 
the control points of fx excludes the origin. It is a well-known property of the 
Bcrnstcin-Bczicr representation of a polynomial / that f{U) lies in the convex 
hull of its control points, where U represents its natural parametric domain. 
Thus, if the hull of the coefficients of a polynomial system excludes the origin, 
this system has no solutions in the parametric domain. We can check whether 
the hull excludes the origin by solving a low-dimensional linear programming 
problem. Megiddo [TS] showed that low-dimensional linear programming prob- 
lems can be solved in linear time (i.e., linear in the number of control points), 
although we have not used Megiddo's algorithm. Note that some other polyno- 
mial bases also have a similar exclusion test; refer to pO] . 

We now proceed to describe our algorithm, the Kantorovich-Test Subdivision 
algorithm or KTS in short. 

Algorithm KTS: 

• Let Q be a queue with [0, 1]"+^ as its only entry. Set S* = 0. 

• Repeat until Q — % 

1. Let X be the hypercube at the front of Q. Remove X from Q. 

2. If X ^ Xe' for all Xe' £ S, 

— Perform the exclusion test on X = B{x^, r) 

— If X fails the exclusion test, 

(a) Perform the Kantorovich test on X 

(b) If X passes the Kantorovich test, 

i. Perform Newton's method on ft,'*'^) , where i is the index 
that satisfies the condition of the Kantorovich test and 
k = xf — starting from to find a zero x* . 

ii. Trace the segment of the solution curve using x* as the 
starting point and going toward x^ + r direction until the 
Xi = x^ + r boundary is reached. 

iii. If the newly found segment is contained in any Xe' € 
S (i.e. the segment has been found before), discard the 
segment. 

iv. Otherwise, compute the new explored region Xe according 
to (HI]). Set S = S\J{Xe}. 

(c) If either X fails the Kantorovich test or X passes the test 
with X % Xe, subdivide X along all n -I- 1 parametric axes 
into 2"+^ equal smaller hypercubes. Add these hypercubes 
to the end of Q. 
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• Check if any two segments of solution curves overlap. If so, remove the 
overlapping part from one of the segments. 

• Join any two segments sharing an endpoint into one continuous curve. 
Repeat until there are no two curves sharing an endpoint. 

A few remarks are needed regarding the description of the KTS algorithm. 

• The subdivision in step 2.c is performed in the case that X passes the Kan- 
torovich test but X % Xe because, in general, passing the Kantorovich 
test does not imply that there is only one solution curve in X. 

• The check in step 2.b.iii is necessary since the segment detected by the 
Kantorovich test may be outside of X. 

• For the same reason as above, certain parts of a solution curve may be 
traced twice and hence must be removed from one of the segments before 
the segments are joined. The overlapping segments can be detected by 
checking if an endpoint of a segment is inside an explored region of another 
segment. The segments sharing an endpoint can also be detected from 
explored regions in a similar manner. Note that there is no ambiguity in 
this step because an explored region, having passed the Kantorovich test, 
cannot contain more than one connected component of the solution. 

• If the Kantorovich test is not applicable for a certain hypercube due to 
the Jacobian of the midpoint being singular, the hypercube is treated as 
if it fails the Kantorovich test and is then subdivided by step 2.c. 

One property of KTS is that it is affinely invariant. In other words, left- 
multiplying / with an n X n matrix A prior to executing KTS does not change 
its behavior. This is the main reason that we introduced icond(/) earlier. 

6 Implementation details 

The implementations of certain steps of KTS are not apparent and thus are 
explained in detail in this section. While this section focuses only on the Bezier 
surface case, all the results herein can be generalized to certain other polynomial 
bases as mentioned in the introduction. 

6.1 Computation of Lipschitz constant 

For simplicity, denote /i^**^-* as h when the choice of (ik) is clear from context. 
The Lipschitz constant for h'{x^)~^h' = g, which is required for the Kantorovich 
test, is obtained from an upper bound over all a; e [—.5, 1.5]"+^ of the derivative 
of g 
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where ^h). (x) denotes the ith entry of [h'{x^) ^h) (x). Let g be the 

Bernstein polynomial that reparametrizes with [0, 1]"+^ the surface defined by 
g over [-.5, 1.5]"+^ (Refer to Section US]) . We have 



max 



g'{x)\\ ^ -. max^ 



a;e[-. 5,1.5]" + ! 2 a;e[Oa]" + i 



2 a;e[Oa]" + i ||j/|| = l 



1 

2 a;eIo7lI" + i' 



n+1 n+1 



^ o ■ _i?ax ^ max jg^^^ (x) \ 



j=i fe=i 



< max max 0,,i,(a;) . 

2 xe[o,i]"+i ^ 

Note that each entry of g' can be written as a Bernstein polynomial efficiently 
because 

= m(Z,_i,„,_i(t) - , (12) 

where Z„i^m_i(t) 

— Zm.m—i{t) — 0, which can be used to compute the control 
points of the derivatives in Bernstein basis from a given Bernstein polynomial 
directly. Hence, the maximum absolute value of the control points of g'^j^. when 
written in Bernstein basis is an upper bound of max^g[o.i]"+i Iffijfel^^)!- Let uj 
denote the Lipschitz constant computed in this manner, that is, 

^ _ (n + 1)2 

o ..,.max ||5jjfe,»i,»2,...,«„+ill' (13) 

ffy7=,^l,^2....,^„+l ^re the control points of 

6.2 The Kantorovich test and solution curve tracing 

Recall that for a hypercube X to pass the Kantorovich test, there must exist 
ani€:{l,2,...,n+l} satisfying ijco < 1/4 and i?(x", p-) C D for all functions 
/j(«fe)'g where fc € — r,x^ + r]. The algorithm, however, need not explicitly 
check the conditions for all values of k. Notice that ui and D are independent of 
k and p_ is an increasing function of rj. For these reasons, KTS only needs to 
check the conditions for the value of k that maximizes 77. Similarly, the explored 
region can be computed solely from the maximizer k. But note also that 77 
is linear in fc, which means that the value of k maximizing tj is either — r or 
x° + r. 

After a hypercube passes the Kantorovich test, the segment of the solution 
curve detected by the test must be traced. Since the Kantorovich test guarantees 
that performing Newton's method on h''^^^ starting on converges for any k € 
[x^ — r, x^ + r], we can trace the segment by repeating Newton's method starting 
on for many different values of k to locate the points on the segment of the 
solution curve. Alternatively, we can perform Newton's method on 
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starting on a;" to find a point on the segment, use as the starting point 
for Newton's method on e > 0, to find the next point on the 

segment, use x"^ as the starting point on ^''+2^) to find the next point on 
the segment, and so on. 

6.3 Reparametrization 

There are two steps of KTS involving reparametrization of polynomials in Bern- 
stein basis, namely the exclusion test and the computation of the Lipschitz con- 
stant for the Kantorovich test. "Reparametrization" in this context means the 
computation of new control points that describe the same function with respect 
to the new parameter domain. Both steps require the reparametrization with 
[0, 1]""*"^ of Bernstein polynomials with n+1 variables, which is a straightforward 
extension of reparametrization with [0, 1]^ of bivariate Bernstein polynomials. 
An example of efficient algorithms for reparametrizing bivariate Bernstein poly- 
nomials with [0,1]^ is discussed in (20]. Alternatively, the polynomials can be 
reparametrized by two applications of the de Casteljau algorithm |16| : one to 
reparametrize the right endpoints with 1 and another to reparametrize the left 
endpoints with 0. 

7 Time complexity analysis 

In this section, we prove a number of theorems leading to the theorem regarding 
the running time of the KTS algorithm. Since both the exclusion test and the 
computation of the Lipschitz constant in the Kantorovich test use the control 
points in their computations, it is useful to find the relationship between the 
control points and the function values of the polynomial defined by them. Recall 
that M is defined as the maximum norm among control points of / and was 
denoted ||/|| earlier. 

We have already shown that f{x) and f'{x) are bounded by AI in ^ and 
([4]). Using the same logic, we can derive a Lipschitz bound on f'{x), i.e., an 
upper bound on /", as follows: 

Wf'ix)- f'iy)\\ <4(n+l)max(mi,...,m„+i)2M|la:-y|l. (14) 

The use of the Kantorovich theorem requires a Lipschitz bound for a slightly 
larger region. If we require a Lipschitz bound for f'{x) over [— e, 1 -t- e]""*"^, 
then we can argue based on the deCasteljau algorithm for evaluating Bezier 
polynomials that 

- f'{y)\\ < 4(n + 1)(1 + e)— (mi,...,m„+0 inax(mi, . . . , m„+i)2M||x - y\\ . 

(15) 

Another useful inequality is that for any x S [0, 1]"^^, 

M||/'(x)t|| > l/(2(n + l)max(mi,...,m„+i)). (16) 
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Equation fTC]) can be established by multiplying both sides of ^ by 
and then using the fact that ||/'(a;)|| • ||/'(a:)t|| > ||/'(a;)/'(a;)t || = ||/|| = 1 for 
any matrix p-norm, and for the infinity-norm in particular. 
Now, we establish a bound that reverses ([3]), namely, 

M<e max ^\\fix)\\ (17) 

for any polynomial /, where 9 is as defined as follows. 

Q ai \ TT I TT max{ mfc - ^ , ^ | 
= 0(mi,...,m„+i) = 11 2^11 T-—7\ 

fc=l \ i=0 i'^i ' ' 



o 



This value of 6* is specific to the choice of the Bernstein-Bezier basis; see [2Uj for 
a discussion of other bases. 

Our proof of (|17p is based on establishing a similar result for univariate 
polynomials as shown by the following lemmas. 

Theorem 7.1 (Srijuntongsiri and Vavasis [lOj)- ^et /(t) he a polynomial system 

fit) = EIlo^'^^™W' o<t<i, 

where bi E R'* . The norm of the coefficients can be bounded by 

ll&^ll <Ci3M^maxJ|/(t)||, 

where 

■^-^ J. J. U — 7 

1=0 j=0,l,...,i-l,i+l,...,m ' 

Lemma 7.2. Let I and h be constants satisfying I < h. Let n be a given positive 
integer. Suppose there exists a function ^(m) satisfying 

||a.|| <eM max ||5(t)|| (18) 

l<t<h 

for any ai € {i = 0,1,..., m) and any univariate polynomial g{t) = 
ai(j)i{t), where <pi(t) denotes the polynomial basis. Suppose also there exists 
a function C{mi, m2, . . . , r7i„) satisfying 

||aii,i2,...,i„ll < C(™l:'7l2, . . . ,TOri) Hiax \\g{xi, X2, . ■ ■ , Xn)\\ (19) 

l<Xi,X2,. ...Xn<h 

for any ai^^i^^,,,^i^ € R'' {ij = 0, 1, . . . ,TOj) and any polynomial in n variables 

mi 7712 rn„ 

g{xi,Xi, ■ • ■ ,a;„) = ^ ^ • • • ^ a^^,^2,...,ir,'P^l{Xl)(t>i2{x2) ■ ■ ■ 0i„(x„). 

il =0 42=0 i„=0 
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Then 



\\bii,...,i„+i\\ < ({mi, 1712, ■ ■ . ,m„)^(m„+i) max ||/(a;i,a;2, . . . ,a;„+i)|| 

l<Xi,X2,---,Xn-\-i<.h 

(20) 

for any G M'' {ij = 0, 1, . . . ,mj), where f is the polynomial in n + 1 

variables defined by 

mi "in + l 

f{xi,X2,---,x„+i)=^--- ^ fei,,...,j„^,(^i,(a;i) • • •0,„+i(x„+i). (21) 

ii— 



Proof. Let f{xi,X2,...,Xn+i) be an arbitrary polynomial defined as in pT 
For any n-tuple /' — {i[,i'2, . . . , i'^) {i'j = 0, 1, . . . , rrij), define 



gi'{xn+i)^ ^ &i;,i^,...,i;,i„+i0i„+i(a;«+i)- 



i„ + l=0 



Applying (fT8|) to g/' yields 



nil < C(™n+i) ,^max ||g/.(a;„+i)|| 



(22) 



for any 6i;,i^,...,j;,j„+i (i„+i = 0, 1, . . . , m„+i). Let Xj, be a point where 
maxi<x„_f_i<h\\gi'ixn+i)\\ is achieved. Define 

kr{xi,X2, • ■ ■ :a;„) = ^ ^ • • • ^ ^ &n,i2,...,i„+i'/'n (a;i)'/'i2 (a;2) ■ • • '/'i„(a;n)0i„+i(a;J')- 

zi— 22— i„— Oz^+i— 

Applying to fc/' yields 



&»i,...,i„+i</>j„+i(a;J/) 



< C(™i,nT-2, ■ • ■ ,"i«) max ||fc//(a;i, . . . ,a;„) 

l^X\ ... .,Xn<h 



(23) 

for any (ii, 12, ... , i„) (ij = 0, 1, . . . , mj). Consequently, by combining (|^^ and 

(Hi, 



< ({mr,+i)\\gr{x*j,)\\ 

C(?7ln+1 



< ^imn+i)Cimi,m2,. ■ ■ ,mn) max ||fc//(xi, . . . ,x„)|| 

^<a:i , . ..,Xn<.h 

< ^(m„+i)C(mi,m2,...,TO„) max ||/(xi, 2:2, . . . , a;„+i)|| . 

l<.Xi,X2,----Xn-\-i<.h 



□ 
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Hence, ([T7)) holds by induction using Theorem 17.11 as the basis and Lemma 
17.21 as the inductive step. 

Recall that the Lipschitz constant uj given by ([T^ is not the smallest Lip- 
schitz constant for h'{x^)~^h over D — [—.5,1.5]""'"^. However, we can show 
that d) < ((n + 1)2/2) Ouj, where u denotes the smallest Lipschitz constant for 
h'{x^)~^h over D. Since uj is computed from the absolute values of the control 
points of g'ijkix), by 

, in + l)\ , , I 

oj < c/max max 0,,t,(a;) 

(n + l)2 I , , M 

= ^ e max max j^^^-fe (a;) | 

< ^" + ^^% max||g'(a;)|| = ^" + (24) 

With this bound on lu, we can now analyze the behavior of the Kantorovich 
test. 

The following is the main technical result of this article. 

Theorem 7.3. Let f{x) — f{xi,X2, ■ ■ ■ ,Xn+i) be a Bernstein polynomial sys- 
tem in n dimensions such that cond(/) < oo. Let be a point in [0, 1]""*"^. Let 
r>0 be such that B{x",r) C [0, 1]"+^ // 

r<(-^)\ (25) 



cond(/) 



whe' 



re 



^° 64(n+ 1)36' • 1.5max(mi,m2,...,m,.+ i) niax(TOi,TO2, . . . ,m„+i)2 ' ''^^^ 

then either 

1. The hypercube B{x^,r) passes the Kantorovich test and the associated ex- 
plored region Xe contains X , or 

2. The hypercube B{x^,r) passes the exclusion test. 

Proof. Let X denote the hypercube B{x^ ,r). Let us introduce additional nota- 
tion for frequently used quantities. Let stand for ||/(a;°)||, and let tt stand for 
||/'(x°)'''||. The proof is divided into two cases by the relative value of the two 
terms in the definition of condition number ^ evaluated at a;°. Let 



Cl 



\J 64(n + 1)36* • (i.5)max(mi,m2,...,m„+i) niax(mi, m2, . . . , TO„+l)2. (27) 



Case 1: ciA/tt < ^jMj^. 
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Note that this case encompasses the possibiUty that = 0, i.e., that x° lies 
on a solution curve. On the other hand, this case requires tt < oo, i.e., 
rank(()/'(a;°)) = n. Note that A/tt is the second term of ([5]) evaluated at 
x° while \/mI(1) is the first term. Therefore, cond{/) > min(M7r, Af/(/i) > 
min(M7r, cfM^TT^). Since ciMtt > 1 and c^Mtt > 1 (combine ^ and it 
follows that the second term dominates the first, hence 

cond(/) > Mil. (28) 

By the hypothesis for this case, cf A/^tt^ < M/cf), i.e., 

< l/(c?Af7r2). (29) 

Let vixP) be the unit-length null vector of f'{x'^), i.e., f {x'^)v{x'^) = 0, 
||w(x°)|| = 1. Let i be such that |tii(a::°)| — 1. Define 



h{x) 



for an arbitrary k S [x^ — r, + r]. By using the facts that ||t;(a;°)| 
|i'i(2^")| = 1, and 



.Vr.^(T)"'^((-^)/V.'. 



v{x^)'^ei 

where Ci denotes the ith column of the identity matrix, it is seen that 

T^^\\h'{x°)-'h{x^\<2{\\nxyf{x^)\\+r), (30) 

for any k £ [x^ — r, x^+r\. (Note that in the infinity norm, \\I ~ vej / {v^ ei)\\ < 2, 
where i is the index such that \\v\\ = \v'^ei\.) The first parenthesized term on 
the right-hand side of ()30|) is clearly bounded by ncf), which in turn is bounded 
by l/(c?Af7r) by 

It follows from ([Ml) , CSl) , and ([261) that co/ cond(/) < 1, hence (co/ cond(/))2 < 
(co/ cond(/)). Combining this with and (pS)) yields 

r <co/cond(/) <co/(Af^). (31) 

Thus, proceeding from ((30)l . 

Note that 4(n l)(1.5)'"^''(™i---'™"+i) max(mi, . . . ,m„+i)2Af is a Lipschitz 
constant of/' on [-5, 1.5]"+i by Consequently, 4(ri-Hl)(1.5)™'''^(™i-'™"+i)- 
max(mi, . . . , 77i„+i)^Af7r is a Lipschitz constant of fix^Yf on [—5,1.5]"+^. 
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From this bound, we derive a Lipschitz constant uj for h'{x^) ^h' over 
[-.5, 1.5]"+i as follows. 

\\h'{x°)-'{h'{y)-h'{z))\\ 
UJ = max 7 7 - 

y,zeD:,y^z WV ~ z\\ 

(i-0^)n.y{f'iy)-nz)) 



y,zGD;y^z \\y - z\ 



2\\f'{xy{riy)-f'{z))\ 

< max 



y,z<ED-y^z lly-^ll 



< 2n ■ max 

y,z£D-y=£z \\y - z\\ 

< 8(n + i)(i.5)max(™i,-,™.+i) max(TOi, . . . ,mn+ifMTT. 



For the last line, we used p5)l . But from (|24|), 

^ < (!i±il!gl^^ < 4(„ + ;L)3g|(^_5)max(mi,...,m„+i) n,ax(mi, . . . , mn+i)^M7r 

(33) 

where w is the Lipschitz constant for h'{x'^) ^h' over [—.5, 1.5]"+^ and uj is the 
computed Lipschitz constant for h'{x^)^^h' over D as defined in 
Combining p2|) with ([33)) yields 



T^cD < 8(71 + i)30(i.5)max(,„i,...,m„+i) jnax(mi, . . . , m„+i)2(co + 1/c?). (34) 

By choice of Cp and ci in and ([77]) respectively, we see that ?/l2) < 1/4 for any 
A: € [x^ — r, a;^ + r] , which is one of the conditions for X to pass the Kantorovich 
test. 

For the other condition, note that \/l — 2h > 1 — 2/i for < ft. < 1/2. 
Therefore, 



1 - VI - 2t]uj 

P- = 

u> 

< 2-q 

< 8(n+l)(co + l/c2)max(TOi,...,m„+i) (by (HID). (35) 

By choice of cq and ci, we conclude that p- < 1/2 and therefore X C 
[— .5, 1.5]"^^^, the domain for which w is a Lipschitz constant. This proves that 
X satisfies the Kantorovich conditions. 
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Finally, the associated explored region Xe contains X because 
1 + VI - 2?7(2; 



U! 



> 



> 



1 



OJ 



1 



4(n + i)35((i,5)max(mi,...,m„+i) niax(mi, . . . , TO„+i)2M7r 
> r, 

for any fee [a;^ — r, + r]. The last line follows from (pS)) and (pij) . 



(by iSD) 



Case 2: ciMtt > ^jMjt^. 

Note that this case encompasses the possibility that tt = cxd, i.e., rank(()/'(a:°)) < 
n. On the other hand, this case requires </> > 0, i.e., is not a root. Since 

cond(/) > min(M7r, Af/(/)) > min(-\/Af/(/)/ci, M/(f), we conclude 



cond(/) > J-2-, 



(36) 



which implies by (|^ that 



,2^2, 



r < 



M 



(37) 



Select an arbitrary x € B(x'^,r). We now derive a bound on /(x) — f{x'^) by 
applying the fundamental theorem of calculus. 



< 



< 



f'{x° +t{x-x°)){x-x°)dt 
\\f'ix°+t{x-x"))\\dt-\\x-x°\\ 



rf 2{n + l)ma.x{mi,...,mn+i)Mdt (by g])) 
Jo 



2r{n + 1) max(mi, . . . , mn+i)M 



< 2{n + 1)cqC^ max(mi, . . . ,m. 



■n+l) 



(by dSZl)) 



< 7^ (by jig) and jl3 



Thus, by definition of i 



Define f{x) such that 



|/(a;)-/(xO)|| < ||/(x")| 



(38) 



/(ii,i2, • ■ • ,ir!+i) = fi2rxi+x1-r,2rx2+X2-r,... 

2rxn+i+x"„^^-r). (39) 
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In other word, / is a polynomial that reparametrizes with [0, 1]"+^ the surface 
defined by / over X. In terms of /, ([55]) is equivalent to 



.f{i)-fii°) 



< 



fix" 



for an arbitrary x e [0, 1]"^^, where x is the rescaled x, and x° is the rescaled 
x" according to ([55]) . In particular, 



max 



fix) -fix") < fix") 



Letgix) = fix)-fix"). By m, 



\c^,,...,^„^,\\ <0- max \\gix) 



for any control point Cj;j^...^i„^j^ of which is equivalent to 



(40) 



(41) 



< 6 ■ max 

*e[o,i]"+^ 



(42) 



for any control point ai-^_,,,^i^_^-^ of fix) (since a constant additive term to a 
polynomial corresponds to a translation of all of its control points) . Substituting 
dH]) into the left-hand side of (gHl) yields 



^21 ,...,2Ti-f 



< 



(43) 



which implies that the convex hull of the control points of fix) does not contain 
the origin. Therefore, X passes the exclusion test. □ 



8 Computational results 

The KTS algorithm is implemented in Matlab and is tested against a number of 
problem instances with three equations and four variables of varying condition 
numbers. Higher dimension problems, especially the ill-conditioned instances, 
require too much computation time due to the large number of hypercubes 
that must be considered. We estimate the condition number by evaluating 
min{l/ 11/(3;") II , II /'(a;^)''' III at the center point x" of every square considered 
by KTS during its execution and also at uniformly sampled points in [0, 1]**. 

Table |3] compares the efficiency of KTS for each test problem with its con- 
dition number. The total number of hypercubes examined by KTS during the 
entire computation, the width of the smallest hypercube among those examined, 
and the maximum number of Newton iterations to converge are reported. Note 
that the high number of Newton iterations of some test cases (the 9th and 10th 
rows of Table [3]) is because the Jacobians of the zeros are ill-conditioned causing 
large roundoff error in the computation of the Newton iterations. 



25 



max{mi, TO2, ms, 7714} 


cond(/) 


Number of 

J.X y v- ui L_/ wo 


Sina,llest 
width 


McLx. Newton 
iterations 


2 


6.60 


641 


.03125 


3 


2 


11.5 


3089 


.01563 


3 


2 


15.5 


673 


.03125 


3 


3 


24.0 


145 


.06250 




3 


50.0 


4273 


.00781 


3 


3 


120 


1009 


.00391 


3 


3 


2.40 X 10^ 


18177 


.00049 


3 


3 


9.88 X W 


15841 


.00195 


4 


3 


1.86 X 10^ 


28881 


.00098 


7 


3 


2.66 X 10^ 


29649 


.00098 


7 



Table 3: Efficiency of KTS algorithm on problems of different condition num- 
bers. 

9 Conclusion and future directions 

We present the KTS algorithm for solving systems of polynomial equations with 
one more unknowns than the number of polynomials. By using the combination 
of subdivision and Kantorovich's theorem, our algorithm can take advantage of 
the quadratic convergence of Newton's method without the problems of diver- 
gence and missing some solutions that commonly occur with Newton's method. 
We also show that the efficiency of KTS has an upper bound that depends 
solely on the condition number of the problem instance. Nevertheless, there are 
a number of questions left unanswered by this article such as 

• Tighter bound on r. Some of the bounds in Section [7] appear loose 
and could potentially underestimate the performance of the algorithm. 
For example, the scalars may be loose, and one step in the argument 
preceding ([3T|) uses the weak bound that x'^ > x since x > 1. Thus, it 
seems like there is room for tightening the analysis. 

A second limitation of our analysis is that we establish a lower bound 
on the smallest hypercube size, which indirectly places an upper bound 
on the total number of hypercubes explored by the KTS algorithm (and 
hence its running time). This upper bound, however, is usually far from 
tight as illustrated by our computational experiments. Thus, a different 
analysis that addresses the number of hypercubes more directly would be 
useful. 

• Using KTS in floating point arithmetic. In the presence of roundoff 
error, we may need to make adjustments for KTS to be able to guarantee 
that the computed intersections are accurate and that all of the solutions 
are found. 
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• Handling singular solutions and degenerate instances. Instances 

containing singular solutions or degeneracy are ill-posed, and our pro- 
posed KTS algorithm does not aim at handling such instances. Certain 
applications, however, look for singular solutions or solutions to degen- 
erate instances. Further investigation on extending KTS to handle these 
situations would be beneficial. 

• Other representations of /. As mentioned in the introduction, we 

assume that / is specified by its (toi + l)(m,2 + 1) • • • (mn+i + 1) Bernstein- 
Bezier control points. In many applications, however, there may be a 
more parsimonious representation. For example, in the SSI problem, two 
surfaces of bi-degree (1711,7112) are separately each represented by (mi -|- 
1)(to2 + 1) control points, hence / = Pi — P2 is fully described by 2(mi -|- 
1) (to2 -|- 1) control points rather than the (mi -|- l)^(m2 -|- 1)^ control points 
needed for the general case. It would be useful if the KTS algorithm could 
work directly on a more concise representation. 

• Extension to general underdetermined polynomial systems. Poly- 
nomial systems with n equations and m > n -I- 1 unknowns generally con- 
tain higher dimension solutions. The subdivision and exclusion test ideas 
still hold for the general case, but a different technique is needed to trace 
an approximation to the intersection surface. 

10 Acknowledgements 

We benefited from a helpful discussion with F. Cucker about condition numbers. 

References 

[1] p. Biirgisser, F. Cucker, and M. Lotz. The probability that a slightly 
perturbed numerical analysis problem is difficult. Math. Comp., 77:1559- 
1583, 2008. 

[2] F. Chaitin-Chatelin and V. Fraysse. Lectures on finite precision computa- 
tions. SIAM, 1996. 

[3] F. Cucker, T. Krick, G. Malajovich, and M. Wschebor. A numerical algo- 
rithm for zero counting. II: Distance to ill-posedness and smoothed analysis. 
To appear, J. Fixed Point Theory Appl, 2009. 

[4] J. Demmel. The probability that a numerical analysis problem is difficult. 
Math. Comp., 50:449-480, 1988. 

[5] P. Deuflhard and G. Heindl. Affine invariant convergence theorems for 
Newton's method and extensions to related methods. SIAM J. Numer. 
Anal., 16:1-10, 1980. 



27 



[6] G. Farin. Curves and Surfaces for CAGD: A Practical Guide. Academic 
Press, 5 edition, 2002. 

[7] R. T. Farouki and T. N. T. Goodman. On the optimal stability of the 
bernstein basis. Mathematics of Computation, 65(216):1553-1566, October 
1996. 

[8] R. T. Farouki and V. T. Rajan. On the numerical condition of polynomials 
in berstein form. Comput. Aided Geom. Des., 4(3):191-216, 1987. 

[9] R. M. Frcund and .J. R. Vera. Condition-based complexity of convex opti- 
mization in conic linear form via the ellipsoid algorithm. SI AM J. Optim., 
10:155-176, 1999. 

[10] G. H. Golub and C. F. Van Loan. Matrix Computations, the Johns Hopkins 
University Press, 3 edition, 1996. 

[11] D. Jiang and N. F. Stewart. Backward error analysis in computational 
geometry. In Marina L. Gavrilova, Osvaldo Gervasi, Vipin Kumar, Chih 
Jeng Kenneth Tan, David Taniar, Antonio Lagana, Youngsong Mun, and 
Hymiseung Choo, editors, ICCSA (1), volume 3980 of Lecture Notes in 
Computer Science, pages 50-59. Springer, 2006. 

[12] L. Kantorovich. On Newton's method for functional equations (Russian). 
Dokl. Akad. Nauk SSSR, 59:1237 1240, 1948. 

[13] P. Koparkar. Surface intersection by switching from recursive subdivision 
to iterative refinement. The Visual Computer, 8:47-63, 1991. 

[14] R. Leary. Global optimization on funneling landscapes. J. Global Opti- 
mization, 18:367-383, 2000. 

[15] N. Megiddo. Linear programming in linear time when the dimension is 
fixed. J. ACM, 31:114-127, 1984. 

[16] N. M. Patrikalakis and T. Maekawa. Shape Interrogation for Computer 
Aided Design and Manufacturing. Springer- Verlag Berlin Heidelberg, Ger- 
many, 2002. 

[17] M. Shub and S. Smale. Complexity of bezout's theorem i: Geometric 
aspects. Journal of the American Mathematical Society, 6(2):459-501, April 
1993. 

[18] D. A. Spielman and S.-H. Teng. Smoothed analysis: an attempt to explain 
the behavior of algorithms in practice. Commun. ACM, 52(10) :76- 84, 2009. 

[19] G. Srijuntongsiri and S. A. Vavasis. A condition number analysis of a line- 
surface intersection algorithm. SIAM Journal on Scientific Computing, 
30(2):1064-1081, 2007. 



28 



[20] G. Srijuntongsiri and S. A. Vavasis. Properties of polynomial bases used 
in a line-surface intersection algorithm. littp://arxiv.org/abs/0707.1515, 
February 2009. 



[21] K.-C. Toh and L. N. Trefethen. Calculation of pseudospectra by the Arnold! 
iteration. SI AM J. Sci. Comput, 17:1-15, 1996. 

[22] D. L. Toth. On ray tracing parametric surfaces. SIGGRAPH Comput. 
Graph., 19(3):171-179, 1985. 

[23] L. N. Trefethen and D. Bau. Numerical linear algebra. SIAM Press, 1997. 



29 



